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Direct numerical simulation (DNS) of turbulent compressible flows is performed using 
a higher-order space-time discontinuous-Galerkin finite-element method. The numerical 
scheme is validated by performing DNS of the evolution of the Taylor-Green vortex and 
turbulent flow in a channel. The higher-order method is shown to provide increased accu- 
racy relative to low-order methods at a given number of degrees of freedom. The turbulent 
flow over a periodic array of hills in a channel is simulated at Reynolds number 10,595 
using an 8th-order scheme in space and a 4th-order scheme in time. These results are 
validated against previous large eddy simulation (LES) results. A preliminary analysis 
provides insight into how these detailed simulations can be used to improve Reynolds- 
averaged Navier-Stokes (RANS) modeling. 

I. Introduction 

Prediction of separated flows about curved bodies remains a significant challenge in computational fluid 
dynamics (CFD). These flows involve large-scale unsteady motion of the separation and reattachment points, 
and are generally poorly predicted by current Reynolds averaged Navier-Stokes (RANS) models. The devel- 
opment of new RANS models for these types of flows will rely on validation with respect to detailed data 
obtained through direct numerical simulation (DNS) on representative test-cases, which capture the main 
flow features. In this work we present DNS results for the separated flow over periodic hills using a new 
unstructured higher-order scheme for compressible high-Reynolds number flows. 

The flow over an array of hills in a channel has been used as a benchmark case for the study of separation. 
It is modeled using a domain which is periodic in the stream-wise and span-wise directions, simplifying the 
boundary conditions in a numerical simulation. Periodic hills were initially examined experimentally by 
Almeida et al. * 1 at Reynolds numbers of up to 60,000 based on the hill height. In this work we consider 
the geometry presented by Mellon et al. 2 This geometry retains the complex flow features from the original 
experiments of Almeida et al. in a smaller domain for the numerical simulation, and the corresponding 
experimental data more closely match the periodic approximation. 

The periodic hills geometry has been investigated both numerically and experimentally over a wide range 
of Reynolds numbers. 2-5 Frohlich et al. 3 performed highly resolved incompressible large eddy simulation 
(LES) at Reynolds number, Re = 10, 595 using two different 2nd-order finite-volume discretizations with 
different subgrid models. Breuer et al. 5 present a comprehensive review of large eddy simulation (LES) and 
DNS performed to date, comparing experimental data with numerical results using DNS up to Re = 5600 
and LES up to Re = 10, 595. The current work extends these studies to DNS at Re = 10, 595. 

Higher-order methods have been shown to be more efficient for simulations requiring high spatial and 
temporal resolution, allowing for solutions with fewer degrees of freedom and lower computational cost than 
traditional second-order CFD methods. 6 As this work is focused on subsonic flows without shocks, the exact 
solution is in C °° , and thus we do not expect the convergence rate of a higher-order scheme to be limited by 
solution irregularity. Spectral methods were employed for the first DNS calculations in periodic domains. 7 
However, as DNS has been performed for increasingly more complex geometries other numerical methods 
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such as finite-difference and finite- volume methods have become more common. In this work, we use a space- 
time discontinuous Galerkin (DG) finite-element method, which extends to an arbitrary order of accuracy. 
Finite-element methods are well suited to problems with complex geometry, where unstructured meshes may 
be required. Higher-order finite-element methods are particularly attractive due to the possibility of using 
local h- and p-adaptation. The use of a space-time formulation allows for adaptation in both the spatial and 
temporal directions, similar to adaptive mesh refinement (AMR) methods with sub-cycling. 

This paper begins with a brief description of the governing equations and our space-time discretization in 
Section II. In Sections III and IV we present initial validation of our numerical scheme by performing direct 
numerical simulation of the Taylor-Green vortex evolution and a fully developed turbulent channel flow. In 
Section V we present results from the direct numerical simulation of the periodic hill problem. Finally, we 
present a summary and discussion in Section VI. 

II. Numerical Method 

The compressible Navier-Stokes equations are written in conservative form as: 


p,t + ( pu z ) }i 

= 0 

( 1 ) 

( pUj),t + ( puiUj +p5ij)'i 

= r b'.i 

( 2 ) 

(pE),t + (pUiH) ti 

= ( TijVj + urTjSij)^ 

(3) 


where p is the density, Uj are the components of the velocity, E is the total energy, p is the static pressure, 
H = E + £ is the total enthalpy, Ty is the viscous stress tensor, kt is the thermal conductivity, T = p/ pR 
is the temperature, and R is the gas constant. The pressure is given by: 

P = (7 - !) ( pE - \pv k v k ) , (4) 

where 7 is the specific heat ratio. The viscous stress tensor, T t j , is given by: 

Tij = p (v i: j + vj'i) - A v k , k Sij (5) 

where p is the viscosity, A = is the bulk viscosity and 5 ij the Kronecker delta. We rewrite (l)-(3) as a 
single vector equation: 

u t + (F I i -FY) !i =0 (6) 

where u = [p, puj , pE] is the conservative state vector, while F{ and FY are, respectively, the inviscid and 
viscous fluxes. Applying a change of variables u — «(v), where v are the entropy variables, the Navier-Stokes 
equations may be rewritten as: 
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where s = log (p/p 1 ) is the entropy. 

We proceed to discretize (7) as follows. The domain, f 2, is partitioned into non-overlapping hexahedral 
elements, k, while the time is partitioned into time intervals (time-slabs), /" = \t n 1 t n+1 ]. Define Vh = 
{u;, w\ K £ [P(k x -0] 5 }, the space-time finite-element space consisting of piece-wise polynomial functions in 
both space and time on each element. We seek a solution v £ Vh which satisfies the weak form: 


E 


| j J.-(wju + W' Xi (F I i -FY)')+ j J wiFjrii - FYrii) + J E(r +1 )«(t” +1 ) - ui(t" )u(F )j = 0 


2 of 20 


American Institute of Aeronautics and Astronautics 


for all w G Vh ■ Here F\rii and m denote numerical flux functions approximating the inviscid and viscous 
fluxes, respectively. In this work, the inviscid flux is computed using the method of Ismail and Roe 9 , while 
the viscous flux is computed using the method of Bassi and Rebay. 10 

The space, Vh, is spanned by the tensor product of ID nodal Lagrange basis functions defined at the 
Gauss-Legendre points. Integrals appearing in (9) are approximated with numerical quadrature using a 
polynomial dealiasing rule with twice as many quadrature points as solution points in each coordinate 
direction. With the choice of basis and quadrature rule, (9) gives a system of nonlinear equations which 
must solved for each time-slab. In this work we use the preconditioned Jacobian- free Newton-Krylov solver 
previously presented. 11 


III. Taylor-Green Vortex Problem 


Initial validation of our numerical method was performed using DNS of the evolution of the Taylor-Green 
vortex. The Taylor-Green vortex evolution is used as a model problem for turbulent flow as it involves only 
periodic boundary conditions, no forcing and a simple initial condition. The flow is solved on an isotropic 
domain, which spans [0, 2-7rL] in each coordinate direction. The initial conditions are given by: 
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where u , v and w are the components of the velocity in the x-, y- and ^-directions, p is the pressure and p 
is the density. The Taylor-Green vortex flow is simulated using the compressible Navier-Stokes equations at 
Mo = 0.1 and Re = 


PqVqL 


= 1600. The flow is initialized to be isothermal (- = — = RTo). 

v p Po u ' 



Figure l. Iso-contours of vorticity magnitude at peak dissipation for the 
Taylor-Green vortex evolution at M = 0.1, Re = 1600, computed 
using 256 3 degrees of freedom. 


Starting from the simple initial condition, the flow becomes turbulent through repeated vortex stretching 
leading to progressively smaller eddies, which are then dissipated to heat through the action of molecular 
viscosity. Figure 1 shows the iso-contours of vorticity magnitude at the point of peak dissipation from a 16th- 
order solution. The figure highlights the presence of small vortical structures captured by the higher-order 
scheme. 

For each simulation the temporal evolution of the kinetic energy 

E k = ^ j \pv-vdVL (14) 
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is monitored. The evolution of the kinetic energy dissipation rate e = — -^r- was computed from the data at 
the space-time quadrature points. Figure 2 shows the temporal evolution of the kinetic energy dissipation 
rate computed with 256 degrees of freedom in each coordinate direction using 2nd- 4th-, 8th-, 16th-order 
spatial discretizations with a 4th-order temporal discretization. The corresponding meshes have, respectively, 
128, 64, 32 and 16 elements in each coordinate direction. Reference data computed from an incompressible 
simulation using a spectral code on a 512 3 grid 12 is also presented. For the 2nd-order scheme, there remains 
significant numerical dissipation at this resolution and the point of peak dissipation is poorly captured. 
With increasing polynomial order, the results relative to the spectral data are significantly improved. The 
dissipation rate for 4th-, 8th- and 16th-order schemes appear to fall directly upon those of the spectral data. 
Zooming in at the point of peak dissipation, Figure 2(b) shows that 8th- and 16th-order schemes do better 
match with the reference spectral data. 
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Figure 2 . Taylor-Green vortex problem at M = 0.1, Re = 1600, computed 
using 256 3 degrees of freedom. 


We present further validation of our numerical method by performing direct numerical simulation (DNS) 
of the Taylor-Green vortex problem at Re = 160. At this lower Reynolds number, we are able to use 
fewer degrees of freedom to resolve the flow so that we may observe the asymptotic convergence rate of our 
numerical scheme. We assess the quality of our numerical solutions by computing individual terms in the 
kinetic energy evolution equation. For compressible flow, the kinetic energy dissipation rate is given by the 
sum of three contributions e = £i + £2 + £3 = — ■ 


e i 

£2 

£3 



2 jj'Sjj s-jj d-} 

I pu k ,kdft 
n 


(15) 

(16) 
(17) 


where Sij = k( u i,j + u j,i) is the strain-rate tensor. Since the flow is nearly incompressible, we expect the 
dissipation due to the bulk viscosity (62) and the pressure-dilatation term (£3) to be small. The kinetic 
energy dissipation rate is then approximately equal to e ss . However, for the compressible simulation this 
does not hold exactly. Time histories of — ^§r, ei, £2 and £3, computed using an 8th-order DG scheme in 
both space and time with 128 degrees of freedom (DOF) in each coordinate direction, are presented in Figure 
2. We note that ei, £2 and £3 must be computed using the “lifted” gradients, accounting for the jumps in the 
solution between elements, in order to be consistent with our DG discretization. Compressibility effects are 
evident in oscillations of the pressure dilatation term (£3), These oscillations are 3 orders of magnitude less 
than the viscous dissipation due to strain (ei), however they do not go away with further mesh refinement, 
but correspond to the true compressible result given the isopycnic initial condition. 
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Figure 3. Kinetic energy dissipation balance for the Taylor-Green vortex 
problem at M = 0.1, Re = 160, computed using N = 8, and 128 3 
degrees of freedom. 
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Figure 4 shows the evolution of ei, e 2 and €3 computed using a the space-time DG scheme with a 

4th-order temporal discretization and 2nd-order spatial discretizations using 48 and 128 degrees of freedom 
(DOF) in each coordinate direction. The pressure-dilatation, € 3 , has a significant bias, contributing large net 
positive kinetic energy dissipation. In particular, at the lower resolution, using 48 DOFs in each direction, the 
pressure-dilatation contribution to the dissipation corresponds to nearly a third of the total. With 128 DOFs 
in each direction, the contribution of the pressure-dilatation term is significantly smaller, while the kinetic 
energy dissipation is predominantly due to the physical dissipation, ei. With increasing mesh refinement, 
the biased pressure-dilatation term decreases toward zero as shown in Figure 5(a). Alternatively, the biased 
pressure-dilatation term is reduced with increasing polynomial order as demonstrated in Figure 5(b). With 
sufficient resolution, the pressure-dilatation term converges to the small (but nonzero) solution presented in 
Figure 3, which exhibits a damped oscillation to zero. The biased contribution of the pressure-dilatation 
term appears to be a numerical artifact of the upwind DG scheme, not present in under-resolved simulations 
using higher-order central-difference schemes . 13 



T T 

(a) N = 2 with 48 3 DOF (b) N = 2 with 128 3 DOF 

Figure 4. Kinetic energy dissipation balance for the Taylor-Green vortex 
problem at M = 0.1, Re = 160, computed using 2nd-order spatial 
discretization. 



(a) N = 2 (b) DOF = 48 


Figure 5. Biased pressure-dilatation term for the Taylor-Green vortex prob- 
lem at M = 0.1, Re = 160, with a) mesh refinement at 2nd-order 
and b) at constant DOF with increasing polynomial order. 
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Our numerical scheme ensures discrete conservation of the state variables: p, put and pE. However, 
secondary conservation laws, corresponding to higher-order moments of the conservative states, such as 
the kinetic energy balance, are not discretely satisfied. The balance of the kinetic energy equation is used 
as a measure of the presence of numerical errors in the simulation without relying on an exact solution. 
Specifically, we compute an error which is the difference between the dissipation integrated over the space- 
time domain and the change in kinetic energy over the given time-interval. The error in the kinetic energy 
as a function of the resolution length-scale h = 1/DOF 1 / 3 is plotted in Figure 6(a). The 2nd-, 4th- and 8th- 
order schemes show convergence of the error at slightly better than formal rate. Increasing the polynomial 
order significantly reduces the error. At the lowest resolution, the 4th-order scheme has an error an order of 
magnitude less than the 2nd-order scheme, while using the 8th-order scheme with the same number of degrees 
of freedom gives nearly another order of magnitude reduction in error. The error versus the corresponding 
CPU time given in terms of TAU benchmark work units 6 in presented in Figure 6(b). The benefit of the 
higher order scheme is clearly seen in this plot. The 4th order scheme on the coarsest mesh is able to attain 
the same error level as the 2nd-order scheme on the finest mesh, but with a cost which is nearly two orders 
of magnitude less. Extrapolating the data, the 8th-order scheme at the lowest resolution requires 4 orders 
of magnitude less work than the 2nd-order scheme to reach a comparable error level. 
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Figure 6. Error in kinetic energy evolution for Taylor-Green vortex problem 

at M = 0.1, Re = 160. 
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IV. Turbulent Channel Flow 

We consider the flow in a channel at Re T = 180, where Re T = is the Reynolds number based on the 
wall shear velocity, u T = \J t w /p, the channel half- width, 5, and the kinematic viscosity, v = p,/ p. This flow 
has been previously studied using incompressible DNS by Kim et al., 14 and provides a good validation of 
our numerical scheme for wall-bounded DNS. The simulations performed are nearly incompressible with a 
Mach number of approximately 0.1 based on the bulk velocity and mean speed of sound. Following Kim et 
al. 14 the domain is of size 47 tS x 26 x 27 t<5, corresponding to stream- wise, normal, and span- wise directions, 
respectively. The domain is periodic in the stream-wise and span-wise directions while adiabatic no-slip 
boundary conditions are applied at the channel walls. The flow is driven by a constant body force applied 
to the stream-wise momentum equation. 

DNS are computed using the space-time DG scheme with a 4th-order temporal discretization and an 
8th-order spatial discretization. A mesh refinement study using spatial discretizations with 96 x 64 x 80, 
144 x 96 x 120, 192 x 128 x 160 and 288 x 192 x 240 degrees of freedom in the stream-wise, wall normal 
and span-wise directions. On the finest mesh, this corresponds to an average spacing in the stream-wise and 
span-wise directions of Ax + « 8 and Az + ss 5 per DOF. The degrees of freedom are clustered towards the 
walls such that the first element in the wall normal direction extends to approximately y + = 5, with average 
spacing in the first element of A y + ss 0.65. The corresponding reference spectral simulation of Kim et al. 
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Figure 7. Channel flow at Re T = 180, instantaneous stream-wise velocity. 

Cross-sections of the full domain of size An5 x 25 x 2ttS. x-, y- and z- 
axes correspond to stream-wise, normal and span-wise directions 
respectively. 


used a mesh with 192 x 128 x 160 DOFs. 

The flow was initialized to the laminar profile on the coarser meshes and turbulence was induced by 
applying the forcing technique of Eswaran and Pope. 15 The forcing is stopped after two eddy turnover 
times, h/u T after which point the turbulence is self-sustaining. The solution was transfered from the coarser 
meshes and used as the initial condition on the finer meshes. Statistics were collected over 20 eddy turnover 
times for all but the 288 x 192 x 240 DOF mesh, for which statistics were collected over 10 eddy turnover 
times. We do not believe that the flow achieved a statistically converged state, however, these simulations 
are sufficient to provide a validation test case for our numerical scheme. 

Figure 7 contains cross-sections of the instantaneous stream- wise velocity on the mesh with 192 x 128 x 160 
degrees of freedom. The cross-sections in the stream- wise and span-wise directions show the presence of large 
eddies on the order of the channel width 28. The presence of streak-like structures near the walls may be seen 
in the cross-section in the wall normal direction. Contours of vorticity magnitude colored by the stream- wise 
velocity are shown in Figure 8. The presence of vortical structures including hairpin vortices are clearly 
visible. 

The mean velocity, u, normalized by u T , using the meshes with 64 and 192 degrees of freedom in the wall 
normal direction as well as the incompressible DNS data of Kim et al. 14 are plotted in Figure 9. Even at 
the lower mesh resolution, the 8th-order DG results show good agreement with the reference DNS results 
obtained using a spectral method. In particular, we note that the viscous sublayer and buffer layers computed 
with the DG scheme match exactly with the DNS data despite this being resolved within only two elements 
in the wall-normal direction on the coarsest mesh. Mismatch of the velocity profile in the core of the flow 
is likely due to not having reached a statistically converged state. Figure 10 shows the corresponding r.m.s. 
velocity fluctuations u rms = ( u'u ) 1 / 2 . The corresponding Reynolds shear stress, uv' is plotted in Figure 
11. The data agrees well with those of Kim et al. with more resolution improving the correlation. 

Again, looking at the balance of higher-order moments of the conservative state provides further validation 
of our numerical scheme. The evolution equation for the Reynolds stresses for incompressible flow are given 
by: 


— — Pij + eij + $ij + Dfj + Dlj + DVj 


(18) 


where denote the material derivative, while the balance terms on the right-hand side of (18) are given 
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Figure 8. Channel flow at Re r = 180, instantaneous contours of vorticity 
magnitude colored by stream-wise velocity. 




(a) 96 x 64 x 80 (b) 288 x 192 x 160 

Figure 9. Channel flow at Re T = 180, mean velocity profile. (Solid lines are 
DG solution, symbols are solution of Kim et al. 14 ) 
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(a) 96 x 64 x 80 


(b) 288 x 192 x 240 


Figure 10. Channel flow at Re r = 180, r.m.s velocity fluctuations. (Solid 
lines are DG solution, symbols are solution of Kim et al. 14 ) 




(a) 96 x 64 x 80 (b) 96 x 64 x 80 




(c) 288 x 192 x 240 (d) 288 x 192 x 240 


Figure ll. Channel flow at Re T = 180, R U v profile. (Solid lines are DG 
solution, symbols are solution of Kim et al. 14 ) 
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We compute the incompressible Reynolds stresses and balance terms in order to compare with previous 
simulations of Kim et al. 14 For this low Mach number flow, the neglected compressible terms are small. 
Figure 12 contains the terms in the evolution of the Reynolds shear stress uv' , normalized by u T , with 
increasing mesh resolution. Overall, the profiles obtained on each mesh are very similar. However, on the 
coarsest mesh large oscillations appear in the pressure-strain and pressure-diffusion terms at the interface 


between elements close to the wall. With increasing mesh resolution these oscillations go away and there 
is excellent agreement with the reference spectral data. We postulate that the large oscillations are due to 
insufficient resolution resulting in non-physical pressure-strain/pressure-diffusion as seen at low resolution in 
the Taylor-Green vortex test case. When present, these numerical artifacts serve as an indication of a lack 
of resolution. 

The statistical distributions of the components of the velocity near the point of peak Reynolds shear stress 
are presented in Figure 13. The stream-wise velocity, u, has an obviously non-Gaussian profile with skewness 
of -0.74 and flatness of 3.50. The wall-normal velocity component, v, deviates somewhat from the Gaussian 
distribution with skewness of 0.28 and flatness of 3.43, while the w has a nearly Gaussian distribution with 
skewness of 0.006 and flatness of 3.26. In the following section we will compare these baseline distributions 
with those seen in different regions of the flow about the periodic hill. 

V. 2D Periodic Hill 

The 2D periodic hill problem has been widely studied as a model problem for separated turbulent flow. 2 " 5 
The geometry we consider was presented by Mellen et al. 2 and has been used as a test case in several studies 
(c.f. Breuer et al. 5 and references therein). The geometry consists of a periodic channel with a 2D hill 
restriction. The size of the domain is 9 h x 3.035 h x 4.5 h in the stream-wise, wall-normal and span-wise 
directions, where h is the height of the hill. The flow is driven by a constant body force to ensure a given 
mass flux. 

We consider the flow at a Reynolds number Re = 10, 595, where the Reynolds number is defined using the 
bulk velocity and the height of the channel above the hill. Figure 14 shows the mean velocity profile, u , and 
the r.m.s. velocity, u rms = {u u') 1 ^ 2 . This flow has several complicated features which pose modeling chal- 
lenges. The flow accelerates up the windward side of the hill, separates over the top and a large recirculation 
region develops on the leeward side of the hill. Reattachment occurs on the flat surface between successive 
hills. Figure 15 shows the instantaneous isocontour corresponding to zero stream-wise velocity colored by 
entropy, which highlights the extent of the reversible flow region. Modeling challenges involve accurately 
predicting the point of separation, which can have significant impact on the size of the recirculation bubble. 5 
Previous LES at Re = 10, 595 have shown the presence of small recirculation regions near the top of the hill 
and at the base of the hill on the windward side. 3 These studies have also noted large span-wise velocity 
fluctuations beyond the reattachment point due to large-scale eddies convecting downstream from the shear 
layer. 3 This phenomenon dictates that significant resolution is needed near the shear layer. 

Simulations were performed at Re = 10,595 using our space-time DG formulation with a 4th-order 
temporal discretization and an 8th-order spatial discretization. We have used a sequence of meshes using 
128 x 64 x 64, 192 x 96 x 96, 256 x 128 x 128 and 384 x 192 x 192 degrees of freedom in stream-wise, wall-normal 
and span-wise directions. At each mesh resolution 8th-order polynomial curvilinear meshes were generated 
to match the solution order by defining the location of the Gauss-Legendre-Lobatto nodes for each element 
using a pseudo-spectral rational interpolation from an underlying structured mesh. The cross-section of the 
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(a) 96 x 64 x 80 (b) 144 x 96 x 120 
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Figure 12 . Channel flow at Re T = 180, balance terms in evolution of 
Reynolds shear stress. (Solid lines are DG solution, symbols 
are solution of Kim et al. 14 ) 
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Figure 13. Channel flow at Re T = 180, distribution of velocity components. 
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(b) llrms 


(a) U 

Figure 14. Mean stream-wise velocity, u and r.m.s. stream-wise velocity, 

Urms = (wV) 1/2 for Periodic Hill at Re = 10,595. 



Figure 15. Isocontour of u = 0 colored by entropy for Periodic Hill at Re = 10, 595. 
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mesh consisting of 16 x 8 x 8 elements, corresponding to the simulation using 128 x 64 x 64 DOFs, is plotted 
in Figure 16. 



Figure 16 . Curvilinear mesh for Periodic Hill. 


The simulations were run for approximately 25 flow through periods. We do not believe this is a suffi- 
ciently long time to achieve a statistical convergence, however, these simulations serve as an initial validation 
of our numerical scheme. As a reference, previous LES simulations of this flow by Frolilich et al. 3 and Breuer 
et al. 5 collected statistics over 55 and 140 flow through times, respectively. 

The mean pressure and wall shear stress on the lower surface of the channel computed on the finest mesh is 
presented in Figure 17. The mean shear stress profile shows the primary separation point at approximately 
x = 0.20, with reattachment at x = 4.37. This matches well with values of x = 0.19 and x = 4.69, 
respectively, for the separation and reattachment points reported by Breuer et al. 5 Table V gives the mean 
separation and reattachment points with increasing mesh resolution. The coarse simulations predict a later 
separation point and earlier reattachment point, under-predicting the size of the separation bubble. With 
mesh refinement the DG simulation results converge to the previously reported data. Detailed analysis of the 
mean shear stress profiles reveals two additional separation points, corresponding to recirculation bubbles 
at the base of the hill on the wind-ward side and at the peak of the hill. The existence of these two smaller 
recirculation bubbles is consistent with previous observations by Breuer et al. 5 




(a) Mean Pressure (b) Mean Wall Shear Stress 

Figure 17 . Mean pressure and wall shear stress on lower surface for the 
periodic hill at Re = 10,595. 


We compare our numerical simulations with LES data of Frolilich et al. 3 Figure 18 contains the 
mean velocity and Reynolds stress profiles computed using the finest mesh resolution at stations with 
x ~ 0.05,2.0,4.0, and 6.0. In the present results, we have plotted the profiles corresponding to data along 
grid lines, as opposed to constant x stations, in order to simplify our post-processing. Figure 18 shows 
good agreement between the present simulations and that of Frolilich et al. Some discrepancies between the 
computed velocity profiles are attributable to a mismatch in the exact location of the plotted data. The 
small difference in the size of the separation bubble computed in our simulation and that of Frohlich et al. 3 
is evident in the velocity profiles at x ~ 4. The Reynolds stress profiles obtained with the DG scheme also 
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Mesh 

Separation Point 

Reattachment Point 

128 x 64 x 64 

0.28 

3.94 

192 x 96 x 96 

0.22 

4.13 

256 x 128 x 128 

0.21 

4.22 

384 x 192 x 192 

0.20 

4.37 

Breuer et al. 5 

0.19 

4.69 


show good agreement with the reference data. However, numerical artifacts are apparent in the profiles, 
which may be due to insufficient statistical convergence. The Reynolds stress profile at x = 0.05 shows large 
values of cross-stream Reynolds stress, w'w'. This “splatting”, corresponding to turbulent kinetic energy 
being transfered from stream-wise and wall-normal fluctuations to span-wise fluctuations, was previously 
described by Frohlich et al. 3 

The corresponding balance terms in the evolution equation for turbulent kinetic energy, k = | u are 
plotted in Figure 19. As with the Reynolds stresses, the computed statistics show good agreement with those 
of Frohlich et al. 3 However, numerical artifacts are apparent in the profiles at the boundary of elements. 
Given similar behavior in under-resolved turbulent channel flow simulations, this suggests that even the 
finest mesh simulation is somewhat under-resolved. 

We compute the probability distribution function of the velocity at five selected points in the domain. The 
location of the selected points are plotted in Figure 20, while corresponding distribution of the components 
of the velocity at each point are presented in Figure 21. Point 1 is selected near the upper surface of the 
channel, away from the hill and is intended as a baseline profile which is similar to the result presented for the 
channel flow. Point 2 is selected over the top of the hill. The presence of “splatting” is observed in the long 
tails of the span- wise velocity distribution relative to the stream- wise and wall-normal velocity distributions. 
At point 3 in the recirculation region, we see mean large negative u and positive v with distributions which 
have tails not significantly longer than the baseline. Points 4 and 5 are in the most energetic part of the 
domain, straddling the mean separation line. At these two points, we see very long tails in the u- and 
u-prohles, corresponding to the large-scale unsteadiness as a result of the separated flow. The skewness of 
both u and v switch sign, showing that they are straddling the shear layer. 

We examine the correlation between the mean strain and Reynolds stresses in order to evaluate the 
validity of the Boussinesq hypothesis: 


Sjj 




‘rru 


(25) 


The contours of the mean strain and the deviatoric part of the Reynolds stresses are given in Figure 22. In 
general, regions of large strain are correlated with regions of large Reynolds stress, however the sign of the 
constant of proportionality in the shear layer changes between the normal and shear stresses. Additionally, 
in the shear layer, the magnitude of the normal strains are similar to the shear strain, unlike an attached 
boundary layer where the magnitude of the normal strains typically is much smaller. Coming over the top of 
the hill, the mean normal strain changes sign while this is not seen in the normal Reynolds stresses. Further, 
at the rise to the hill, a large normal stress is present in a region with no corresponding strain. Hence, a linear 
eddy viscosity model will over-predict strongly the turbulent diffusion in this region. These observations are 
indicative of the insight that can be gained from these simulations to improve RANS modeling. 


VI. Summary and Conclusions 

In this paper we have performed direct numerical simulations of turbulent compressible flows at low mach 
number using a higher-order space-time discontinuous-Galerkin method. We validated our numerical scheme 
for turbulent flow by studying the Taylor-Green vortex problem. We demonstrated that the upwind DG 
scheme gives a net positive bias in the pressure-strain contribution to the kinetic energy dissipation for under- 
resolved simulations. This biased contribution disappears with increased resolution and we demonstrated 
formal error convergence rates up to 8th-order. 

Next, we validated our numerical scheme for wall-bounded turbulent flows by performing DNS of turbu- 
lent channel flow at Re T = 180. We presented results for mean flow and higher moment quantities which 
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Velocity Profile 


(a) Velocity Profile, x = 0.05 



Reynolds Stresses 

(b) Reynolds Stress, x = 0.05 



Velocity Profile 

(c) Velocity Profile, x = 2.0 



Velocity Profile 

(e) Velocity Profile, x = 4.0 



Reynolds Stresses 

(d) Reynolds Stress, x = 2.0 



Reynolds Stresses 

(f) Reynolds Stress, x = 4.0 




Velocity Profile Reynolds Stresses 


(g) Velocity Profile, x = 6.0 (h) Reynolds Stress, x = 6.0 

Figure 18. Velocity and Reynolds stress profiles at various stations in the 
domain for the periodic hill at Re = 10,595. (Solid lines are DG 
solution, symbols are solution of Frohlich et al. 3 ) 


16 of 20 


American Institute of Aeronautics and Astronautics 


0.2 


0.1 

0 


- 0.1 


- 0.2 



Dissipation 

Production 

Pressure-Strain 
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(a) X = 0.05 


(b) x = 2.0 




(c) x = 4.0 (d) x = 6.0 

Figure 19. Balance terms in evolution equation for turbulent kinetic energy 
at various stations in the domain for the periodic hill at Re = 
10, 595. (Solid lines are DG solution, symbols are solution of 
Frohlich et al. 3 ) 



Figure 20. Location of points for computing velocity distributions on the 
periodic hill at Re = 10, 595 
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(a) u at point 1: (—0.35,3.29) (b) v at point 1: (—0.24,4.81) 


(c) w at point 1: (0.07,3.73) 



(d) u at point 2: (—0.12,4.23) 


(e) v at point 2: (—0.18,4.01) 


(f) w at point 2: (—0.03, 3.00) 



(g) u at point 3: (0.25,3.70) 


(h) v at point 3: (—0.43, 5.72) 


(i) w at point 3: (—0.13,3.5) 



(j) u at point 4: (0.33, 3.00) 


(k) v at point 4: (—0.38,3.31) 


(1) w at point 4: (0.01,3.40) 



-1 0 1 2 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 


U V w 

(m) u at point 5: (—0.31,2.79) (n) v at point 5: (0.23,2.84) (o) w at point 5: (—0.02,3.21) 

Figure 21 . Velocity distributions at selected points for the periodic hill at 
Re = 10, 595. (Values in parenthesis correspond to skewness and 
flatness). 
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Figure 22 . Reynolds stress and mean strain contours for the periodic hill 
at Re = 10, 595. (Contour lines correspond 30 uniform intervals 
between [-2,2] for strain and [-0.04,0.04] for Reynolds stresses.) 
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were consistent with results previously presented in the literature. 

Finally, we performed direct numerical simulation of turbulent flow in a channel with periodic hill con- 
strictions. We validated our results by comparing with previously computed LES simulations. We have 
presented detailed statistics showing velocity distributions at various regions of the flow. Analysis of the 
resulting flow fields could provide insights that would facilitate improvements to RANS models. Future work 
will involve performing higher-fidelity simulations over larger times to provide statistically converged data, 
data. 
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